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ABSTRACT 



(N ■ 
> ■ 

' Context. FUV radiation strongly affects the physical and chemical state of molecular clouds, from protoplanetary disks to entire galaxies. 
Aims. The solution of the FUV radiative transfer equation can be complicated if the most relevant radiative processes such us dust scattering 
and gas line absorption are included, and have realistic (non-uniform) properties, i.e. if optical properties are depth dependent. 
Methods. We have extended the spherical harmonics method to solve for the FUV radiation field in externally or internally illuminated clouds 
taking into account gas absorption and coherent, nonconservative and anisotropic scattering by dust grains. The new formulation has been 
implemented in the Meudon PDR code and thus it will be publicly available. 

Results. Our formalism allows us to consistently include: (;) varying dust populations and (ii) gas lines in the FUV radiative transfer. The FUV 
^ ^, penetration depth rises for increasing dust albedo and anisotropy of the scattered radiation (e.g. when grains grow towards cloud interiors). 
I Conclusions. Illustrative models of illuminated clouds where only the dust populations are varied confirm earlier predictions for the FUV 
O penetration in diffuse clouds (Aj/<1). For denser and more embedded sources (Av>l) we show that the FUV radiation field inside the cloud can 
differ by orders of magnitude depending on the grain properties and growth. Our models reveal significant differences regarding the resulting 
^ physical and chemical structures for steep vs. flat extinction curves towards molecular clouds. In particular, we show that the photochemical 
. . , and thermal gradients can be very different depending on grain growth. Therefore, the assumption of uniform dust properties and averaged 
^ ' extinction curves can be a crude approximation to determine the resulting scattering properties, prevailing chemistry and atomic/molecular 
abundances in ISM clouds or protoplanetary disks. 
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1. Introduction 



Far-UV (FUV) radiation (hv <13.6 eV) strongly affects the physical and chemical state of dusty molecular clouds in many 
evolutionary stages: from star forming regions ([Lequeux et al. 198T] IStutzki et al. 19881 |Bally et al. 1998[) and protoplanetary 



disks dJohnstone et al. 19981 lAikawa et al. 2002ll. to circumstellar envelopes around evolved stars (Huggins & Glassgold 1982 



Habing 1996 1 and supernova remnants (Shull & McKee 1979NChevalier & Fransson 19941 ). Thus, the accurate knowledge of the 



intensity of the FUV radiation field as a function of cloud depth is of crucial importance in a plethora of astrophysical environ- 
ments. Penetration of FUV radiation strongly depends on dust grains properties through the scattering of photons, but it also 
depends on the gas properties (chemical composition, density, etc.) through the absorption of hundreds of discret electronic lines 
from the most abundant species (H, H2, and CO). This proccess is, in addition, an efficient excitation mechanism for molecu- 
lar species (B lack & van Dishoeck 19871 Sternberg & Dalgarno 19891. Gas absorption lines reach extremely large opacities and, 
due to saturation, they can be very broad and fully absorb the FUV continuum. 

The so called spherical harmonics method, in which the specific intensity of the FUV radiation field is expanded into se- 
ries of Legendre polynomials, is an efficient way to solve the plane-parallel radiative transfer equation if gas opacity is ne- 
glected and if dust grains have uniform optical properties, e.g. the same extinction cross-section, albedo and scattering phase 
function ( Flannery et al. 1980 Roberge 1983[ l. Nevertheless, astronomical observations over the full spectral domain show a 
more complex scenario, where dust grain populations evolve depending on the environmental conditions from polycyclic aro- 
matic hydrocarbons (PAHs) and very small grains (VSGs) to bigger grains (BGs) likely formed by accretion or coagulation 
([Boulanger et al. 1988||Desert et al. 19901 iJoblin et al. 19921 IDraine 20031 IDartois 200511. Also, the average extinction law {e.g.. 
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|Cardelli, Clayton & Mathis 1989| l is based on observations toward low-extinction line of sights (Ay < 5), and it has been ques- 
tioned by recent observations toward more embedded regions (Ay > 15). A better knowledge of the extinction properties at large 
Ay is critical. In particular, there is evidence that the reddening curve tends to flatten at high extinction depths (iMoore et al. 2005T l. 
consistent with grain growth and dust processing along the line of sight. Therefore, the attenuation of FUV radiation will dra- 
matically depend on the (generally poorly understood) grain composition and optical properties that, of course, are likely to 
change from source to source according to the interstellar (ISM) and circumstellar (CSM) dust life-cycle. In addition to this 
dust-shielding, self-shielding through gas line absorption can result in an efficient protection of H2 and CO, and the starting 
point of a rich chemistry even in irradiated media such as protoplanetary disks, translucent clouds, starbursts galaxies or, more 
generally, photodissociation regions (PDRs; see lHollenbach & Tielens 19971 for a review). 

The spherical harmonics method also has been implemented to study the radiative transfer and dust extinction in galax- 
ies as a whole by associating the source function with the emissivity of a given distribution of stars through the galaxy 
jdi Bartolomeo et al. 19951 |Baes & Dejonghe 2001| ). Uniform grain properties and the absence of gas line absorption are as- 
sumed. For unidimensional problems, the spherical harmonics method is found to be by far the most efficient way to solve for 
the radiative transfer equation compared to Monte Carlo or ray tracing techniques (Baes & Dejonghe 2Q 0j|. 



The detailed information provided by high angular resolution observations (e.g.l Gerin et al. 2005. iGoicoechea et al. 20061 1. 
revealing fine differences even between similar sources, should be followed by a sophistication in the radiative transfer modeling. 
Inclusion of gas (discrete line absorption) and varying grain populations (e.g. different extinction curves) as a function of cloud 



depth requires a modification of the original method (Flannery et al. 1980 Roberge 1983 1. In this work we present an extension 



of the spherical harmonics method for a radiative transfer equation with depth dependent coefficients in plane-parallel geometry. 
We used this method to solve for the radiation field in illuminated clouds at wavelengths longer than Lyman cut-off at ~912 A 
taking into account gas absorption and scattering by dust grains. The method can also include the source function for embedded 
emission of photons, and therefore it can explicitly take into account any source of internal radiation. 

In Sees. 2 and 3 we present the formulation of the method while in Sees. 4 and 5 we show several astrophysical applications 
to understand the role of FUV penetration for the photochemistry of molecular clouds. In particular, we present a few examples 
including H Lyman lines, H2 electronic transitions within the Lyman and Werner bands and CO electronic transitions together 
with varying dust properties. The penetration of FUV radiation for the typical conditions prevailing in a diffuse cloud (such us 
^ Ophiuchi) and in higher extinction objects (such as the Orion Bar or a strongly illuminated protoplanetary disk) are discussed. 

2. The equation of radiative transfer with variable coefficients 

The specific intensity of radiation, /^(i,ju), in plane-parallel geometry is a solution of the radiative transfer equation: 
dh{s,n) cr,(s) , , , 

= -[a'i(i) H- o-,i(i)]/,)(i,//) H- — — R,x{s,^i,iJ.)h{s,n)dn + j^is) (1) 

OS 2 J_i 

where the spatial scale s and the angle 9 - cos^^jj are the independent variables and where the dependence of quantities on 
wavelength A and on s has been explicitly written. In the most general problem, a^is) = a^^is) + a^(.s) is the line-plus-continuum 
absorption coefficient, o-^{s) is the dust scattering coefficient, jA{s) is the emission coefficient of any source of internal radiation 
and iJ^(i,/i,/i') is the angular redistribution function (we assume that the radiation field has azimuthal symmetry about normal 
rays). In this work, the opacity is due to coherent (no energy redistribution in the scattered photons), nonconservative (a fraction 
of photons are absorbed), anisotropic scattering by dust grains as well as to gas line absorption, that is: 

dr - -(a^ + cr^) ds (2) 
(note that t increases toward the decreasing direction of s) and the radiative transfer Eq. ([T]i gets the more familiar form 



dh(T,iJ.) to,i(T) r+' , , , . 

j" ^Ia(t,p) — j RA{T,IJ.,p)lA(T,p)dfl - S A(T,fi) ^ lA(T,fi) - Sa(t,Ii) 



(3) 



where oja - 7^77^ is a new effective albedo (the dust scattering cross-section over the total dust-Hgas extinction cross-section) 
which tends to the pure dust albedo for wavelengths free of lines, but tends to (true gas absorption) at the line cores. Intermediate 
values are found in the line wings. S*^ - ^ ^'^^^ is the source function for the true emission by "embedded photon sources". In 
the following we assume that 5^=0. Thus we ignore dust thermal emission (negligible in the FUV for ISM clouds) or any other 
source of internal illumination. Hence, our source function only corresponds to the external illumination photons scattered by 
dust grains. However, inclusion of 5^ in our method is trivial. The interested reader is refererred to Appendix lAl 

The cloud extends from t = to t = T,„nv with a possibility that T^ax - Boundary conditions require Ia(j,ij) to match the 
incident intensity at t = and t - T„,ax- Note the implicit sign convention on fi: 9 - n points towards positive values of r, that 
is yu = -1 for a ray perpendicular to the cloud and penetrating into it from t = (see Fig. [TJ. Thus, boundary conditions specify 
X'in) - /i(0,//) (jU < 0) wdx^ifA - hijmax^lJ) ifi > 0), where ;^'='(yu) are the illuminating radiation fields reaching both cloud 
surfaces (of course they can be different). 
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Fig. 1. Adopted geometry and sign conventions for a cloud with embedded sources of photons 5*(yu) and illuminated at both 
surfaces x^^P)- 



Compared to other works where the spherical harmonics method has been applied to solve for the FUV radiation field (e.g. 
Flannery et al. 1980, ^Roberge 1983[ Idi Bartolomeo et al. 1995. ,Baes & Dejonghe 2001] ILe Petit et al. 2006l l. the optical prop- 



erties in the radiative transfer equation (e.g., effective albedo and asymmetry parameter) are wavelength- and cloud depth- 
dependent for the first time. 

3. The spherical harmonics method for line and continuum transfer 

3A. The Pl approximation 

In this method, the angular dependence of the radiation field /(t, yu) is expanded in a truncated series of Legendre polynomials 
Piifi) which form a complete orthogonal set within the range (-1,1) in which fi varies: 

L 

/(t,//) = ^(2Z+1)/,(t)P,(^) (4) 

1=0 

where the dependence on A is no longer shown. In the following sections we show that the mean intensity of the radiation field at 
each depth point 7(r) has the simple form J(t) = /o(t), i.e. the first coefficient of the expansion in Eq. (HJl, which is often the only 
quantity needed for the integration of radiation field-dependent physical parameters (e.g. photoionization and photodissociation 
rates). This is one of the reasons why the method is so attractive. However, a large number of expansion terms has to be used 
in order to correctly sample the angular dependence of the radiation field, we typically use L + I = 2M - 20 (note that dust 
scattering can be highly anisotropic at the considered wavelengths). 

If the grain scattering phase function p(t, cos 0) only depends on the angle between the incident and scattered radiation, 
R(T,fj.,fi') can also be expanded (see e.g.. lChandrasekhar 19601 Roberge 1983 1 as: 



R(T,fi,fi') = ^(2/+ l)cr,(T)P,(//)P,(//') (5) 

/=o 

in terms of the cr/(T) coefficients of the Legendre expansion of p(t, cos 0): 

L 

P(t, cos 0) = Yj{21 + 1) cr,(T) P,{cos 0) (6) 



The standard model of scattering by interstellar grains ( Henyey & Greenstein 1941 1 assumes the simple scattering phase function: 



1-^' 

picas®) = ^ — —3^ (7) 

(1 + g - 2g cos&yi^ 

which can be also expanded in Legendre polynomials in terms of the "g-asymmetry parameter" (=< cos& >) i.e., the mean 

r+l 

angle of the scattered radiation (g - 1 /2 J ^ jj. pip) d/j., with p - cos ©). Here we adopt a r-dependent Henyey-Greenstein phase 
function (other phase functions can be used if they can also be expanded). Therefore we write: 

L 

piT,cos@) = ^(2Z+ l)g,iT)P,icos®) (8) 



where giij) - criir) and goir) - 1. Thus, the angular redistribution function RiT,p,fi') has two obvious limiting cases, gij) - 
(isotropic scattering) and gir) = +1 with RiT,p,fi') = dip + p') (pure backward or forward scattering). 
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Substitution of Eqs. (|4]i and Q into the transfer equation (O and using appropriate recurrence formulae leads to the finite 
(L + 1) set of coupled, linear, first order differential equations in the unknown //(t) coefficients, with I - 0, L. 



IfUr) + DfM = {21 +!)[!- cn{T)^i{T)] f,(T) 



(9) 



where /' - df/dr. We recall that compared to Roberge (1983) this is not a constant coefficient equation so numerical integration 
is necessary. In the "Pi approximation" a sufficiently large odcf^ L value has to be chosen to obtain an accurate solution of the 
problem. The system (|9]l can be written as: 



f'(T)=A-l(T)f(T) 

with: 





(10) 



A(t) 





2/!-' 



2/1-1 



3/1-1 



{L-l)h-^_, 




(11) 



where: 

hl{T)^(2l+\)(\-Cj(T)cr,(T)) 



(12) 



In summary, we have to solve for a linear boundary value problem with non constant coefficients with the additional difficulty 
of huge variations of the total opacit}3 within small variations in the wavelength and cloud position grids, e.g. from A in a. 
saturated line center (t^ ~10^) to A in an adjacent (line free) continuum region (r^ ~10). In the following, we show an extension 
of the spectral method of Flannery et al. (1980) and Roberge (1983) to solve for the FUV radiative transfer. 



3.2. The eigenvalues solution 
3.2.1. Numerical solution 

The A-1(t) matrix has L + 1 = 2M eigenvalues which are real, non-zero and non-degenerate and which occur in positive-negative 
pairs, see Appendix A of Roberge (1983). Using a similar notation as Roberge, let kmir), m - ±1, ■ ■ ■ , +M be the eigenvalues 
verifying k^mij) - -Ki{t), and R(t) be the matrix of eigenvectors, that is: 

A-\t),jRj,„(t) = kUT)RlmiT) (13) 

j 

which also verifies the Ri-m(j) - i-^)' Rim{T) relation. The depth-dependence of the eigenvalues kmif) and eigenvectors Rimir) 
complicates the solution of the problem compared to the (only dust) problem with uniform optical grain properties. The compu- 
tation of A;„,(t) and Ri,„iT) is given in AppendixICl 

The R(r) matrix of eigenvectors can still be used to define an auxiliary set of variables y(r) - R-'(T)f(T), or: 

-1 M 

/,(t) = ^ Rl,n(T)y,n{T) + J] Rhnir) JnAr) (14) 

-M 1 

SO that 

f'=A-iRy (15) 
Therefore, in terms of the new y(r) variables, Eq. (fTOb can be rewritten as: 

y' = R 'A 'Ry-R iR'y = Ky-R iR'y (16) 

To write Eq. ( fT6l ) we have used the fact that (R-^A^^R)/,,, - ki di,„ and thus K(r) is a diagonal matrix with the k,„{T) eigenvalues 
of A-*(t) on its diagonal. The fact that R'(t) adds the last matrix term in Eq. ( fT6b due to the depth-dependence of the 



' For even values of L, A is singular (e.g. Roberge 1983) 

^ We also developed the formalism to solve Eq. jlOi through ^n;/e differences JAscher et al. 1995) . For only dust continuum transfer, results 
are almost identical (within ~0. 1%) to those obtained with the spherical harmonics method (which is ~2 times faster). However, when line 
absorption is included, the finite difference numerical solution always oscillate at the core of saturated lines and no optimal solution is found. 



Goicoechea & Le Bourlot: The penetration of FUV radiation 



5 



coefficients. This term is neglected in Le Petit et al. (2006). However, R' is not null neither when the grain optical properties 
depend on the cloud depth (even if gas is neglected) nor when gas line absorption is included (even if grain properties are 
uniform). Unfortunately, the system of Eqs. ( fT6l ) is uncoupled only if the R' y term is null (as in Roberge 1983), otherwise 
more manipulations are required to solve the problem consistently. If we define Q - R 'R'y = -Ly, then Eq. (fTSl l can be 
simply written as: 



y'm = l^m(T)y,n + [L y]„,(T) 



(17) 



for m - ±1, +M. In order to solve this particular problem we turn the system of differential equations ( [TT] ) into an integral 
problem. To do that we first introduce the following integral equation: 



Cjfi "t" 



y]m(f) dt 



(18) 



where am(T) is an arbitrary function so that am{Tm)-Q. The system of Eqs. ( fTSl l represents a general set of integral equations that 
verify ym(Tm) - Cm (to be found from the boundary conditions). If a given function ym is a solution of the above equation, by 
taking its derivative with respect to r one gets: 



Jm W = a'm(T)ym(T) + [L y]„,(T) 



(19) 



which means that y^ as defined in Eqs. dlSl l is also a solution of the original system of differential Eqs. ( fTTI i if and only if 
a'^ir) - km(T). Therefore, fl,„(r) - k,„(t) dt. This demonstration shows that the k,„ eigenvalues of A"' (and no others) are the 
right exponential factors that do attenuate the radiation field, which is consistent with the original problem described by Eqs. ( fTOl l. 
In the present work we solve Eqs. ( fTST i with an iterative scheme^ and thus compute: 



r -/' 



k„(t')dr' 



lLf">lnit)dt 



(20) 



by using an appropriate (physical) initial guess for yj,"', where n is the iteration step. This iterative procedure shows that the 
solution if forced, at any step, by the exponential factor e^,,, ''«•(''> follow the behavior dictated by the "true" eigenvalues of the 
problem (i.e. those of the original coupling matrix A"') that are known before the iteration procedure is started. In Appendix IB] 
we give details on the error bound associated with the iterative scheme and we show that the numerical solution derived for the 
FUV radiation field correctly satisfies the original system of Eqs. ( fTOl i. 

At each iteration step we have to compute the integration constants C,„ by a convenient selection of t,„. To ensure that only 
exponentials with negative arguments appear, it is necessary to set r„, = for m < and t,„ - T,„ax for m > 0. In order to have 
easier to read equations, we now introduce some convenient notations: 



E^{t) = exp k,n(t) dt\ (m < 0) 



or 



E,„(t) = exp 



(m > 0) 



Note that Ezjj) = E;„(t), and E;„(0) = 1. We also define: 



E^{t) = exp 



(m < 0) 



or 



= exp 



(m > 0) 



with ^^(Tmfl.t) - 1 and E^j(t) x EJt) = £,1,(0) - E^(Tmax)- Using the above notations, we have: 

qm(t) dt (m < 0) 

E^(t) 



y,„iT) = E^,{t)C,„ + 



Jr E;„it) 



q,„(t)dt (m>0) 



(21) 



(22) 



(23) 



(24) 



Note the change of sign in the second equation due to the inversion of . To further simplify these expressions, we define: 

^m(T)= r^f^qjt)dt (m<0) (25) 
Jo EJt) 

K(r)^ f"'°^^^q.„(t)dt (m>0) (26) 
Jt E + {t) 



which satisfy D^(Q) - and D^(T,„ax) = 0. Therefore, the y(T) variables are finally written compactly as: 
y„,(T) = £™(t) C„, - D-Jt) (m < 0) 



(27) 
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y„,(T) = C„ + DliT) (m > 0) (28) 

and the original //(t) terms in the Legendre expansion of the radiation field /(t, ji) are then given by: 

-1 M 

fl(T) = Yj ('^'« '^-(^^ - '^m(^)) + Z ^'"'^^^ '^-^^^ + '^'i'^^^) ■ ^^^^ 

m=-M m=l 

3.2.2. Boundary conditions: Clouds witli two sides illumination 

We consider a unidimensional plane-parallel cloud of finite size with an external radiation field at both cloud surfaces (t = and 
T = Tmax) defined hy x^(l-i) and;^'^(yu) respectively (see Fig. [Til. From Eq. (|29] l we have: 

-1 M 

fi(0) = 2 R,,„(Q) C,n + Yj ^'."(0) (Cn, + £>:,(0)) (30) 

m=—M m=l 

-1 M 

/Kw) = 2 ^„„(w) (C,„£Jw)-Z);(w)) + 2^-(w)C.. (31) 

m--M m— 1 

At the T = side, the solution must match, at each A, the incoming radiation field with yu < 0, i.e. 1(0, /S) - x'(m), with 

L 

7(0, /z) = ^(2Z + l)/,(0) W or : (32) 

/=o 

-I L M L 

7(0,//) = 2 ^(2/+ l)7?aO)7^;(yu) + 2(C,„£,:(0) + D;(0)) ^(2/ + l)/;„„(O)P,0«). (33) 

m=-M 1=0 m=l /=0 

At the T = T„,fli side, the solution must match, at each A, the incoming radiation field with yu > 0, i.e. I{t„u,x7P) - X^itA^ with: 

L 

7(w,//) = ^(2/ + 1) /,( w) P,(jt/) or : (34) 

/=() 

-I CO ML 

)-7)„-(T„,,,)) _^(2/ + 1)7?/,„(t„,„) PKa*) + ^ C„, ^(2/ + 1) 7?,,„( w) P,(//). (35) 

m=-M /=() m=l /=() 

Nevertheless, since the order L of the expansions is finite, the boundary conditions 1(0, i-i) = x (t^) ™d I(T„ax,fj) - X^(t^) can 
not be satisfied at all /j angles. In this work we use Marcfe'^JI conditions that require 7(0, yU < 0) and 7(r„,ai,ju > 0) to match the 
incident radiation fields at L + 1 = 2M strategic angles //, given by PL+\(Mi) - 0, that is, the roots of the Legendre polynomial of 
degree L+ 1. Note that in these (/ = ±1, +M) angles, the solution of the radiation field I(T,iUi) is "exact". 
To further simplify the boundary conditions relations, we now define: 

L 

Ti,„(0,iud = ^(2/ + 1) 7;ta(0) P,(jjd (fii < 0) (36) 

/=o 

L 

Tim(Tmax,l^i) = ^(2/ + i) Rlm(rmax) Pl(l^i) (Mi > 0) (37) 

/=0 

which gives: 

-1 M 



1(0, Hi) = <^'« ^-^0, fid + ^ (C„, (0) + D„+ (0)) r,„,(0, fid (38) 

m=-M m=l 

-1 M 
^(Tmax, f^i) — ^ ^ ^ iCmE^(Tfj^x) 7)^^(t^^ ^-)) T im(Tmax, f^i) "t" ^ ^ ^ im(Tmux, f^i) (39) 



m=l 



See e.g., Sen & Wilson (1990) for a different choice of boundary conditions. 
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Therefore, the desired constants at each iteration step are solutions of the 2M x 2M linear system (m = excluded): 



M 



2 BunC,„^Hi (40) 
m - -M 

with the B„„ coefficients as define in Table [l] and where 



riO,iJd - 2:,1f=i £>;(0) r,M(0,A*,) < O) 

I^(T„,ax, fii) + Ii,7,L-M D„(T,„ax) Ti,„(T„jax, Mi) (Mi > 0) 





m < 


m > 


yU, < 


7'™(0,yU,) 




yU, > 







Table 1. B„„ coefficients for the two sides illumination boundary conditions in Eq. (l40l i 



For semi-infinite clouds (Tmax-°°) with only one side illumination at t=0 (mi < 0), boundary conditions have to be modified 
to take into account the no radiation condition at t-co {jj^i > 0). It is straightforward to show that the C,„ constants are then 
solutions of the same linear system shown in Eq. (l40l i with the Bi^ coefficients now defined as in Table|2]and: 

(1(0, Mi) (Mi<0) 

"'-[O (Mi>0) ^^^^ 



Bill, - 


m < 


m > 


A', <0 


T.iMfii) 





Mi > 





T imi.^ max 1 f^i) 



Table 2. coefficients for the one side illumination boundary conditions in Eq. (|40] | 



3.3. Iterative procedure 

At very large optical depths (e.g. deep inside the cloud or at the core of saturated lines) the intensity of the radiation field tends 
to zero. Hence, the simplest way to initiate the iterative process is to set Q - R 'R' y = 0. However, this may be far from the 
real solution, and more realistic guesses should be tried. In practice, the assumption r — > oo may be too crude and one can add 
the effect of the external radiation perpendicular to the cloud that penetrates deepest in the cloud, i.e. attenuated by the smallest 
eigenvalue k-t\ (that associated with the radiation field in the \m\ - 1 direction). Thus, we guess a first set of ymij), that we call 
y%{T), from the linear system: 

-1 M 

m=-M m=l 

with 

= i 7(0,-1) exp[-^i(T) r] + ^ 7(t,„«,, 1) exp[^_i(r)(T,„„ - t)] (44) 

Note that only the / - terms have to be considered. As noted by Flannery et al. (1980) and Roberge (1983), the presence of dust 
scattering implies that 4^ 1, i.e. the radiation field attenuation factor at large depths is not simply e but dominated by 
the e ^''^ factor This conclusion obviously applies for the present case with the difference that k±i is now depth-dependent and 
includes line absorption. This important result can modify the intensity of the FUV radiation field inside optically thick clouds by 
orders of magnitude depending on the dust grain optical properties. At lower optical depths (e.g., diffuse clouds), the attenuation 
factor still contains an important contribution from additional terms (A;±2, ^±3, ■■■)■ 

Now that we have an educated guess for the ym{T) variables, we can estimate the new term in Eq. ( fT6] l carrying the depth- 
dependence of the gas and dust coefficients, i.e. the Q = R"' R y term. Note that R^'R' need to be evaluated only once, so 
numerical cost is limited. However, special care should be taken for the R' derivation. Details of the R"^ inversion and R' 
derivation are given in AppendixiDl 
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We briefly describe the iterative computation of Q: we start by using Q* = R 'R' Y" and then compute a first set of C", from 
the boundary conditions. With these first C,*^ and Q" variables we can now use the general expression Eq. ( |20] i to compute a new 
set of ymir) to derive a more refined Q term, and start this proccess again until some prescribed level of convergence in Q is 
reached. Thus, if « is the iteration index, Q("+'* is computed from Q("+'> - R^' R' Y^"^'^ with: 



y] 



(T) = Cri^£„(T)-D;(«)(T) (m<0) 



Those expressions have to be computed at each iteration by numerical integration. 



(45) 
(46) 



3.4. Mean intensity and FUV photon escape probability 



Once we have obtained the full depth and angular description of the intensity of the radiation field /(t,jU/) through the /;(t) 
coefficients, we show here the simple form that J(t) takes. The angular average of the specific intensity is defined as: 



J(r) 



1 r+i 

^ 2 j 1 ^^'^''^^^^ 



From the expansion of I(t, fi) we have: 

X+I /-+! 



(47) 



(48) 



J Po(//) dji - 2. Therefore, as anticipated in Sec. 13.11 the mean intensity of the 
radiation field at each wavelength and depth reduces to J(t) - faij), that is: 



(49) 



where we use the fact that Rq,„(t) - 1 for all m and t. Despite the simplicity of this relation, in many cases of astrophysical 
interest (e.g. a two sides illuminated cloud) one needs to distinguish the fraction of the radiation field coming from each side of 
the cloud. In this case, two half sums have to be computed. In Appendix|E]we give the analytic formulae to compute the mean 
radiation intensity /"^(t) coming from each side. The resulting J'^ir) values can be used to evaluate the escape probably of any 
FUV photon emitted within the cloud, e.g. within H2 line cascades. In particular, the probability for a photon emitted at t-t' 
(inside the cloud) to reach t=0 (or T=Tmax) is given by the 7"(r')/7"(0) (or J^{t')I J^ijmax)) intensity ratios. These probabilities 
can then be further used to determine the H2 level detailed balance. We also note that in this method the first terms of the intensity 
expansion in Eq. (01 are directly related to the moments of the radiation field, i.e. /o(t) - Jij), the mean intensity; /i(t) - H(t), 
the Eddington flux; and /2(t) - 3 K{t) - J(t) where K(t) is the A'-moment. 

From the numerical point of view, the methodology described in the previous sections has been implemented in the 
Meudon PDR cod^ a photochemical model of a unidimensional plane-parallel stationary PDR ( ILe Bourlot et al. 1993t 
ILe Petit et al. 20061 and references therein) and will be the FUV radiative transfer method used in the code. In the following 
sections we illustrate several of the new possibilities with some relevant astrophysical examples. 




1050 iim nso 

vfRvclcnpLii [a] 



Fig. 2. Radiative transfer models for a cloud with a total extinction of Ay—l and a density of nn-lO cm , illuminated at both 
sides by the mean ISRF. Part of the resulting FUV spectra (~912-1300 A) close to the cloud surface is shown. The blue spectra 
correspond to a model with R' =0 in Eq. (fT6l ) (depth dependence neglected), and the red one corresponds to the new "exact" 
computation. 



* Available at http://aristote.obspm.fr/MIS/ 
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4. Applications: Comparison with previous approaches 

In this section we compare the main differences of the new exact computation versus the Hne-plus-continuum approach (R' - 0) 
used by Le Petit et al. (2006) in the Meudon PDR code. Since the previous version of the code used a single-dust albedo and 
g-asymmetry parameter with no wavelength or depth dependence, and the extinction curve was not related to the grain properties 
used in the model, here we just make the comparison by assuming R' = in the new computation, and limit ourselves to the 
uniform dust properties case. In the following examples we explicitely include all the H, H2 and CO electronic absorption lines 
arising from rotational levels up to J -6 (for H2) and 7=1 (for CO). The FGK approximation (IFederman et al. 19791) is applied for 
the rest of levels. Note that the exact method allows one to take into account the overlaps between H, H2 and CO lines neglected 
in more crude approaches. 

Apart from having a radiative transfer method to consistently solve for the dust grain varying populations problem (Section|5]l, 
the next largest difference between the new computation compared to Le Petit et al. (2006) is the effect of line-wing absorption 
of back-scattered radiation. At line core wavelengths, photons penetrating into the cloud are purely absorbed by the gas (the 
effective albedo equals 0). Due to saturation and opacity broadening, many absorption hnes become very wide deep inside 
the cloud. As a consequence, the FUV radiation field is more attenuated than in the (only) dust continuum transfer case. At 
continuum wavelengths free of lines, a fraction of photons coming from the external illumination sources can be absorbed by 
the dust (depending on the exact dust albedo value) or be back-scattered (depending on the exact g value) and provide an 
additional contribution to the radiation field at the cloud surface (about 10% of increase for g -< cosO >^ 0.6). At line wing 
wavelengths, where dust and gas opacities are of the same order (and the effective albedo is in between and the grain albedo), 
some of the back-scattered photons can again reach the surface of the cloud while another fraction will be absorbed in the wings. 
Therefore, as shown by our calculations, line wings are "numerically more challenging". The fraction of absorbed photons in 
the line wings depends on the wavelength separation to the line core and on the transition upper level life time (because it 
determines the resulting line profile broadening). To illustrate these differences we consider a cloud with a constant density 
nH=10^ cm"^ and a total extinction depth of Av=l, illuminated at both sides by the mean interstellar radiation field (ISRF, 
X - 1) as defined by Draine (1978)| These physical conditions resemble those of a diffuse cloud such as parts of f Ophiuchi (e.g. 



Black & Dalgarno 1977 1. An uniform grain size distribution similar to that of Mathis et al. (1977) is assumed. Figure|2]shows part 
of the resulting FUV spectra (~9 12- 1300 A) close to the cloud surface. These spectra clearly show that the effect of H2 line wing 
absorption of back-scattered photons is larger in the exact computation compared to the R' = approach jLe Petit et al. 2006] l. 
Note that this is true only for H2 lines. Atomic hydrogen lines exhibit the opposite effect, i.e. a decrease of the line wing absorption 
of back-scattered photons compared to the R' = approach. Figure [TO] shows the impact of the same two, exact and R' - 0, 
computations in the resulting cloud structure (left: H/H2 transition and right: H2 photodissociation rate). In spite of the different 
line profiles predicted by each type of model, the final cloud physical conditions remain very similar. Therefore, we conclude 
that all computations made with the previous version of the Meudon PDR code jLe Petit et al. 200 6*). where line transfer was 
computed (assuming R' = and uniform dust properties), are consistent with the present exact calculation. The larger effect of 
the H2 line-wing absorptions in the exact calculation increases the attenuation of the illuminating radiation field, which results 
in a H/H2 transition layer slightly shifted to lower extinction depths. This general result obviously applies to any FUV radiative 
transfer model including gas line absorption compared to (only dust) continuum models, i.e. the contribution of gas absorption 
(H2 lines mostly) decrease the photoionization rate (of neutral carbon particularly) and the photodissociation rate of species 
with thresholds close to the Lyman cut. An adventage of including gas line absorption is that predicted spectra can be directly 
compared with spectral observations provided by FUV telescopes. 




10 * 10 ^ 10 ^ 10 ^ 0.1 0.2 0.3 0.4 0.5 

Ay Ay 



Fig. 3. Impact of the new 'exact' radiative transfer computation compared to an alternative approach that assumes R'-O 
jLe Petit et al. 2006b . Grain properties are uniform in all the cloud (MRN). Left panel: H/H2 transition. Right panel: H2 pho- 
todissociation rate as a function of cloud depth. A cloud with a density of nn^lQ^ cnT^, a total extinction depth of Ai/=1 and 
illuminated at both sides by the mean ISRF is considered. These results show that for the case of uniform dust grain properties 
the error associated with R'=0 assumption is small. 
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5. Applications: Grain growth, varying dust populations 

With the method presented in Sec.[3l we can now consistently explore the effect of more realistic (non-uniform) dust properties 
in the FUV penetration into more embedded objects e.g., dense molecular clouds or protoplanetary disks. As a representative 
example, we present several models for a dense and strongly illuminated cloud (with an ionization parameter of y Inn - 1 cm^) 
with grain radii varying dust populations. From the chemical point of view we only concentrate here on the effects that the 
different FUV attenuation depths have on the classical H/H2 and C^/C/CO layered structures predicted by PDR models. In 
particular, we consider a cloud with a constant density hh - n{H) + 2 n{H2)-lQ^ cmT^ and a total extinction of Av-2Q which is 
illuminated at both sides by 10^' times the ISRF. These physical conditions resemble those of a dense PDR such as the Orion Bar 
(e.g. lTielens & Hollenbach 1985l l or a photoevaporating disk around a massive star (e.g. Johnstone, Hollenbach, & Bally 1998). 
At any depth we consider that dust grains follow a size distribution dn - da given by: 

riair) = ^ riajir) = ^ A,(r) nuiT) a''^' da «/ -(r) < fl(r) < fl,>(T) (50) 

where a± refers to the grain radius distribution lower and upper limits and / = 1, ...,n refers to each component of the grain 
mixture. In Eq. dSOl l we have explicitly particularized for the simple power-law case, although more complicated problems may 
require other prescriptions of ng (e.g. such as those in 'Weingartner & Draine 2001 1. Grain properties were taken from Laor & 



Draine (1993) for silicates and graphite. With these tabulations we compute the optical parameters of the grain mixture for each 
wavelength and cloud depth. In particular, we compute the Qats, Qsca and Qe^i efficiencies and the grain albedo QscalQe.xt- We 
finally use an g^-asymmetry factor averaged over the grain distribution as (see e.g.. lWolfire & Cassinelli 19861) : 

Y^iJ^cp- gi(a. A, t) Qscaia, A, r) «a,;(T) da 

^-^w = -^"^^^ : (51) 

X ^ fl QscJa, A, r) ria/j) da 

Afterwards, the extinction curve A{A, t)IAy{t) and the absolute dust extinction coefficient a^^if) are determined at each depth 
and used to settle the total line-plus-continuum opacity (as defined in Eq. |2]) and the effective albedo. The dust extinction 
coefficient (cm"') is given by ad{A, t) - ng na^ Qexi, where is the number of dust grains (per cm-'). Thus, we compute: 



na^ V e^,,(fl,T)A,(T)nH(T)fl^ 



Q\^,{a,T)Ai{T)nH{T)al'\da (52) 



The A,(r) grain coefficients are determined at each depth position assuming that the gas-to-dust mass ratio has to be constant 
(~100) in the whole cloud (i.e., the number of grains is reduced if grain sizes increase). However, in order to keep the grain 




Fig. 4. Grains mixture optical parameters as a function of wavelength and cloud depth (the red curve corresponds to the illumi- 
nated cloud edge Ay=0 and the blue curve to the center of the cloud at Ky-IQ) for the "MRNto BGs" (left) and "VSGs toMRN" 
(right) examples respectively. The shaded region shows the spectral region taken into account in the FUV radiative transfer 
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mixture homogeneous, the Asu/^Cra ratio is kept fixed. Contribution of discrete absorption lines, i.e. the contribution of a^Cr), 
is included in similar fashion as described in Le Petit et al. (2006; Sec. 4.3). The total opacity at each depth is then given by: 



a«\A(A) 



t^T,i = 1 1 + I — ; — dTv = 



^^^l\dry (53) 



where all the variables are depth dependent and where we have assumed that, in the visible band, the extinction is only produced 
by dust and therefore we use dT,i - dry - ds to relate the spatial scale with extinction depth. Note that we compute the 
extinction curve, at each cloud position, directly from the derived grain properties. 

For this "grain growth example" we consider that grain radii increase as a function of the cloud depth according to: 

aijT) = GiJO) + [aijTc) - fl;,±(0)] (^j (54) 

where a,,±(0) defines the grain radii at the edge of the cloud (t = 0) and a,,±(Tc) refers to the grain radii at the center of the cloud. 
We chose ^±-2/3. Obviously, this is just an illustrative example since we do not explicitly solve for the grain nucleation/growth 
(e.g. Salpeter 1974) nor the erosion/sputtering problem (e.g. lBarlow 19781 1. which depends on the particular type of source. The 



crucial point here is to provide a method to consistently solve for the FUV radiative equation if, as suggested by observations, 
the grains size distribution changes toward embedded objects jMoore et al. 20051 ) and/or if spatial fluctuations of the gas to dust 
ratio do exist along the line of sight dPadoan et al. 2006] l. 

In the following, grains follow a power-law distribution of sizes given by Eq. ( l50l ) with yS,=3.5 at each cloud position. A 
mixture of silicates and graphite grains defined by Asu/Acra = 1-1 and with fl_(Av - 0)-5 nm, a+(Av - 0)=250 nm and 
a-(Av - 10)=50 nm, a+(Av = 10)=2500 nm was selected. Therefore, the grain mixture at Av=0 corresponds to the size 
distribution proposed by Mathis, Rumpl & Nordsieck (1977; this size distribution is called MRN hereafter) to fit the mean 



galactic extinction curve (see also Fitzpatrick & Massa 1990 1. At Ay^lO grains have grown by a factor 10 and we call them 
Big Grains (BGs). In the second example we only change the size distribution to fl_(Av - 0)-l nm, a+(Av = 0)=50 nm and 
a-(Av = 10)=5 nm, a+(Ay = 10)=250 nm. Thus, the grain mixture at Av=0 corresponds to very small grains (VSGs). At Av=10 
grains have grown by a factor of 5 and follow a MRN distribution again. The third final example considers a uniform grain size 
distribution (MRN) in the whole cloud. The resulting optical properties, extinction curves, dust opacities. A, coefficients and radii 
distributions for these examples are shown in Figs.|4l|5]|6]and|2](fe/f panel), respectively. 

Some time ago, Sandell & Mattila (1975) emphasized that the albedo and anisotropy of dust grain scattering have important 
effects on photodissociation rates for ISM molecules. The present computation of the FUV radiation field (continuum+lines) 
at each cloud position (see Fig. [8] for the resulting FUV spectra at different Ay) allows an explicit integration of consistent C 
photoionization rates together with H2 and CO photodissociation rates. Once the FUV radiation field has been determined and 
the photo rates calculated, steady-state chemical abundances are computed for a given network of chemical reactions. Finally, we 
compute the thermal structure of the cloud by solving the balance between the most important gas heating and cooling processes 
dLe Bourlot et al. 1993ilLe Petit et al. 2006l and references therein). 

Depending on the grain properties these examples show FUV radiation fields that change by orders of magnitude at large 
Ay (Fig. |7] right panel). Note that the mean radiation intensity at the cloud surface 7(0) cannot be larger than the illuminat- 
ing radiation field itself, i.e., J(Q)/x < 1- The exact ratio depends on the particular dust scattering properties (~0. 53-0. 54 for 
these models of optically thick clouds). The influence of the different grain distributions in the attenuation of FUV radiation 
is obvious, the FUV penetration depth is larger when dust scattering is more efficient, i.e., when grain albedo and scattering 
anisotropy increase (as dust grains grow toward bigger grains). Note that the only difference between models is the change of 
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Fig. 6. Dust mass absorption coefficients (per gas gr) as a function of wavelength and cloud depth (the red curve corresponds 
to the cloud edge Av=0 and the blue curve to Ay-lO) for the "MRN to BGs" (left) and "VSGs to MRN" (right) examples 
respectively. The different grain material A, coefficients required to keep a constant gas-to-dust mass ratio are also shown as a 
function of Ay in the small insets. 
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Fig. 7. Left: Adopted grain averaged radii distribution for the "MRN to BGs" (left) and "VSGs to MRN" (right) examples 
respectively. Right: Resulting mean intensity of the FUV continuum (at ~ 11 32 A) as a function of the cloud depth for the three 
different varying grain populations discussed in the text. The ordinate shows the mean intensity normalized by the illuminating 
radiation field (y=10^ in Draine's units). 



grain size distributions across the cloud. Therefore, the assumption of uniform dust properties and averaged extinction curves 
can be one of the crudest approximations made to determine the resulting cloud physical and chemical state. Figure |9] shows 
the impact of the different grain growth curves on the resulting cloud structure: kinetic temperature, H2 photodissociation rate, 
C photoionization rate and CO photodissociation rates (left column), H/H2 transition, and C^/C/CO abundances (right column). 
The different intensities of the FUV radiation field for each dust population result in very different photoionization and photodis- 
sociation rates which ultimately determine the prevailing chemistry. This conclusion qualitatively agrees with earlier calculations 
for ISM diffuse clouds (Roberge, Da lgarno & Flannery 1981) l and should be extended to more embedded objects where there are 
observational evidences (e.g. lMoore et al. 2005l l of flatter extinction curves (consistent with grain growth). The H/H2 and C^/C 
layered structures in our models are different even in similar sources (same density and illumination) if grain properties signif- 
icantly disagree, or if dust grains vary along the observed region. Different ionization fractions, molecular ions enhancements, 
and C^/C/CO abundances should thus be observed. In particular, photochemistry can still be important at large Ay if anisotropic 
scattering of the illuminating radiation is efficient (e.g., "MRN to BGs" model). In this case, CO photodissociation and carbon 
ionization still dominate the CO destruction and C^ formation respectively deeper inside the cloud. As a result, the predicted 
abundance of neutral and ionized carbon at Av=10 is enhanced compared to standard MRN dust models (see Fig.|9]l. 

Secondly, the intensity of the FUV radiation field also determines much of the thermal structure of the cloud through the 
efficiency of the grain photoelectric effect, the dominant heating mechanism (e.g. Draine 1978). Since FUV radiation penetrates 
deepest when dust grains are bigger, the photoelectric heating rate is kept high deeper inside the cloud. Thus, a larger fraction 
of the gas is maintained warm at large extinction depths. Warmer temperatures also affect the rates of chemical reactions with 
activation energy barriers. For the smallest dust grains, FUV attenuation is so high that photoelectric heating soon becomes 
inefficient and the gas is colder at large extinctions depths. Note that since grain ionization is very large in the surface of the 
cloud (due to the high illumination in the selected example), the maximum efficiency of the photoelectric effect, i.e. the maximum 
temperature, is reached deeper inside the cloud where the grain ionization has decreased. The general effects described here must 
play a significant role in illuminated sources where grain growth takes place, specially in protoplanetary disks, circumstellar 
envelopes around evolved stars and dense molecular clouds near Hn regions. In these cases, the FUV penetration depth is 
increased if dust grains evolve toward bigger grains, leading to larger photochemically active regions. 
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Fig. 8. Radiative transfer models for a cloud with a total extinction of Av=20 and a density of nn-iO^ cm"^, illuminated at both 
sides by 10^ times the mean ISRF. Part of the resulting FUV spectra (~912-1300 A) at different extinction depths: Ay=0 (cloud 
surface), Av-O-l, Av=2 and Av=3 are shown in each box. In each panel, the red (blue) curve corresponds to the "MRN to BGs" 
{"VSGs to MRN") example. 



Conversely, molecular species such as CO will be more abundant in irradiated regions where the smallest grains dominate the 
extinction efficiency. FigurefTOlshows the effects of grain growth in a diffuse cloud (Ai/=1), with a density of hh-IQ^ cmT^, and 
illuminated by the mean ISRF. Although the resulting variations are not so large compared to optically thick clouds, the different 
photoionization and photodissociation rates also translate into different atomic and molecular abundances. 

In particular, the C^/C and C^/CO abundance ratios change up to a factor ~10 depending on the assumed grain properties. 
Note that for optically thin clouds, the mean intensity at one surface can have a significant contribution from the other side 
illumination (that increases with the scattering efficiency). As an example, the mean intensity at Ay=0 in the "MRN to BGs" 
grain model {J{Q)lx-Q-6'i; red curves in Fig.fTOli is a factor ~20% larger than in the "VSGs to MRN" model (blue curves). This 
effect slightly modifies the dissociation and ionization rates at the cloud surface. 

In summary, as gas photodissociation and heating determine much of the chemistry in FUV irradiated gas, the resulting 
source structure is severely altered by the assumed (or observed) grain properties. Therefore, understanding dust properties and 
grain variations in individual sources is a crucial step to determine the source physical and chemical state. 
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Fig. 9. Impact of the different FUV radiative transfer models on the kinetic temperature, H2 photodissociation rate, C photoion- 
ization rate and CO photodissociation rate (left column), H/H2 transition and C^/C/CO column densities (right column). A cloud 
with a density of nH=10^ cm"^, a total extinction depth of Av-2Q and illuminated at both sides by 10^ times the mean ISRF is 
considered. Although not clearly seen in these boxes, all physical parameters show an horizontal tangent at Av = 10, consistent 
with their null variation with respect to the depth position at half cloud (as expected for a symmetrically illuminated cloud). 
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Fig. 10. Same as Figure|9]for a cloud with riH-lO^ cm ^, a total extinction depth of Ay^l and illuminated at both sides by the 
the mean ISRF. Dust grains grow according to Eq. ( |54| ) with A^= 1.086 Ty=0.5. 
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6. Summary and conclusions 

1 . An extension of the spherical harmonics method to solve for the radiative transfer equation with depth dependent coefficients 
in plane-parallel geometry has been presented. The method can be used to solve for the FUV radiation field in externally or 
internally illuminated clouds, taking into account gas absorption and coherent, nonconservative and anisotropic scattering by 
dust grains. Our extended formulation thus allows to consistently include (/) gas fines and (ii) varying dust populations. 

2. We have shown that the penetration of FUV radiation is heavily influenced by dust properties. According to the dust ISM 
and CSM life-cycle, such properties likely change from source to source but also they change within the same object. The 
FUV penetration depth rises for increasing dust albedo and anisotropy of the scattered radiation when grains grow at large 
Ay (as suggested observationally). Therefore, the modeled physical and chemical state of illuminated molecular clouds, 
protoplanetary disks or entire galaxies can be altered by large factors if a more reafistic treatment of the interaction between 
radiation and matter is considered. 

3. The new formulation has been implemented in the Meudon PDR code and thus it will be publicly available. Particular 
examples where only the dust populations are changed show intensities of the FUV radiation field that diff'er by orders 
of magnitude at large Ay. Therefore, the resulting photochemical and thermal structures of molecular clouds can be very 
difl'erent depending on the assumed grain properties and growth. 
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Appendix A: Inclusion of embedded sources of emission {S* i^O) 

In this appendix we give the recipe to include the true emission by "embedded sources of photons" in the method described in 
section|3] In this case the source function includes the scattering of photons by dust grains plus a non null S \ - term (see 
Eq.[3]and Fig.[Tli, where y^(s) is the emission coefficient (line or continuum) of any source of internal radiation. 

Firstly, the angular dependence of S *(r, fj) has to be also expanded in a truncated series of Legendre polynomials Pi(^) as; 

L 

S*{T,fi) = J]^2l+l)si(T)P,(M) (A.l) 

1=0 

where the dependence with A is omitted. The inclusion of Eqs. ( lA.ll ) into the transfer equation (O leads to an additional term in 
the set of coupled, linear, first order differential equations in the unknown //(r) coefficients: 

V/_i(r) + (/ + l)/,;i(r) = (21 + 1) [1 - oj(T)cr,(T)] /,(t) - (21 + 1) s,(t) (A.2) 
Therefore, the system of equations ( IA.2I ) is now non-homogeneous and can be written as: 

f'(T) = A-\t) f(T) + A-\t) g(T) (A.3) 

where gi(T) = -si(t)/(1 - w(t) cr;(r)). Although the method can be easily used for anisotropic source functions, in most practical 
applications, the embedded sources of photons emit isotropically and therefore the terms in the expansion of S*(t) in Eq. ( lA.ll i 
reduce to s/(t) - S*(t)6k), and thus, ^/(t) reduces to go(T) = -5'*(r)/(l - w(r)) with ^;(t) = if / 0. Using the same set of 
auxiliary variables y(T) - R"'(T)f(r), Eq. ( fTSI ) is now written as: 

y' = Ky-R-'R'y + KR-ig (A.4) 

Note that by inserting RR"' between A"' and g, we have simplified R"' A"' g as KR"' g. This result is particularly usefujf], 
since it avoids computing A ' completely. Hence, the last matrix product, Q = KR ' g, makes the system non-homogeneous: 

y'„ = k„(T) y„ + [L y]„,(T) + qm(T) (A.5) 

Eq. dA.Sb can also be solved with the iterative scheme described in section[3]by including the additional term, i.e.. 



^ ill 



" + e-L ^'-C')'"' ([Ly(«^],„(f) + q,n(t))dt 



(A.6) 



It is straightforward to show that the //(t) terms in the Legendre expansion of the radiation field I(t, fi) are still given by Eq. 
The only change compared to the 5'*=0 case is that the qm(T) variables in the D'^(t) and D^„(t) integrals (Eqs.l25]andl26ll have to 
be substituted by ^m(T) - 'qmi^), that is: 

Jo -fcmiO 

{ ""^ iq,n(t) -q,n(t)) dt (m > 0) (A.8) 

The iterative procedure can now be initiated taking into account that at large optical depths the intensity of the radiation field is 
isotropic and tends to the ratio of the true emission to the true absorption: 

I(T ^ cx,) ^ (A.9) 
1 - u(t) 

In practice, the assumption r — > oo may be too crude. We have computed that by adding the effect of the external radiation that 
penetrates deepest into the cloud, the iterative scheme is more robust. Therefore, the first set of ym(T) variables in the iterative 
procedure, y^„(j), are computed from the linear system: 

+M S (t) \ 1 

y Rim(T)yl(T) = f^PP'^^T) = ^"^7 + -7(0,-1) exp[-fci(T) t] + -/(t™„ 1) exp[^_i(T)(w - t)] (A.IO) 

where only the / = terms are considered. 

We have successfully applied the above method by associating 5 * to thermal emission of dust. These kind of computations 
are useful if the radiative transfer calculation is extended to the IR domain (A > \ fim), where scattering of IR photons by dust 
grains is still significative. In the FUV domain, S * can represent any source of internal illumination. In a future paper we plan to 
include "secondary" line photons in the embedded source function. This line FUV radiation field arises from the H2 radiative de- 
excitations that follow the H2 excitation by collisions with electrons and cosmic rays dPrasad & Tarafdar 19831 ) and is generally 
poorly treated. However, it constributes to molecular photodissociation deep inside molecular clouds where the continuum FUV 
radiation field has been attenuated. 



This is true whatever the isotropy properties of the source function are, and not only for the isotropic case. 
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Appendix B: Numerical solution and error limits 

In Sec. 13.2.11 we turn the system of differential equations (fTTl ) into and integral problem (Eqs. [TSl l that we solve numerically 
through an iterative scheme (Eqs. l20t . In this appendix we provide a bound on the error associated with this procedure and we 
verify that the derived solution satisfies the original system of equations (fTOl i. 

(n) 

Given a numerical approximation y;,, to the true solution y,„, we investigate if our iterative proccess converges for all A and 
Ay of the wavelength and cloud depth grids. Thus, we compute: 

with [Ly(«)]^^ (f) = Lmiit) yf\t) (B. 1) 

and write the error in step n + 1 as A^^'' - - y™^'' where m- ±1, +M. Note that this is the difference between the true ym 
(unknown) and our numerical approximation at step « + 1. Since the above equations are linear, aJ"^'' reduces to: 



>,(«+!) 



(t) 



c\ 



(n+l) 



-a,„(t) 1 



Aif i'(T) - e' 



-a„,(t) T 



j(t)Af\t)dt\ 



(B.2) 



because the boundary conditions term (C,„ - Cj;"^'*) is small and damped almost everywhere by the exponential term (as shown 
numerically). If we now define /^^"^•^'^^ - max, | A'"'(r)|, the maximum error at iteration step n in the Legendre expansion of order 
! (/ = 0, 1, L) at any depth position, then: 



+M 



By taking the maximum error at iteration step n at any depth position and at any Legendre order, A*"''*''*^ - max,,, AJ^'"''^'*^, we 
arrive to a severe upper limit to the error between the true solution and the numerical approximation at step n+l: 

^(„+l),MAZ < ^(n).MAX I T [«„(r)-«,„W] j = /^(n),MAX . ^ 4) 

i=-M ^-''"». ' 

Therefore, convergence is guaranteed if ^ < 1 as A'^"^'* *^'^^ < A*'' '*^'*^ Obviously convergence occurs also for less restrictive 
conditions but this is harder to constrain. In our computations we find Jl<\ for almost all wavelengths and depth positions. Only 
at some specific locations in the {Ay) grid (those associated with some line wings), J?l can take values < 10. However, a close 
look at successive variations of A^^'^ - aJ^' at those locations shows that A*,"^'' - A*^' is effectively null after a few iterations. 

A final test to validate our numerical solution is to compute the numerical derivative of our solution and compare f with A f ' 
(see Eq. (fTOb ). Although grain properties are kept uniform, inclusion of gas absorption makes R'(t) and thus L'(t) 
in Eq. ( fTTI i. Figure IbTI shows a typical example for a test cloud with Ay = 1 and uh - 300 cm"^, illuminated by the standard 
radiation field on both sides. In particular, we compare the /; (/=0) component of f with At' alA- 914.26 A, a H2 line wing with 
a total optical depth of 80. Hence, variations of physical conditions along the cloud are large. It can be seen that the agreement is 
excellent. In a continuum "free of lines" wavelength range, agreement is perfect, and there is nothing to show. Hence, the derived 
numerical solution is a very good approximated solution to the radiative transfer problem. 



Appendix C: Eigenvalues and eigenvectors of A ^t) 

We describe here our method to compute the eigenvalues and eigenvectors of the A"*(r) matrix (see Eq.(fT3]l). Note that A(r) and 
A"'(t) have the same eigenvectors, but A;,7/(t) and k,„{T) eigenvalues respectively. 

A first trick is to turn this diagonalization problem to a symmetric problem. Let us call R„,(t) an eigenvector of A(t) with 
Rimij) components and A;,7,'(t) eigenvalues. Thus, R(r) is the matrix formed by the R,„(r) eigenvectors and we can write: 







2/1- ' 



2/22' 



3/13I 



{L-\)h-^\ 




' Ron, ' 




Ri)m 


Rim 




Rim 


Rim 




Rim 








RL-l,m 




RL-l,m 


■. Rlih , 




K Run 1 
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with hiir) — {21 + 1) (1 - oj{t)o-i{t)). If we now define G(t) as diagonal matrix with ^//(t) - hj^^ir), left-muhiphcation of the 
previous equation by G(t) and insertion of the identity matrix I = G"'(t) G(t) between A(t) and R(t) give^: 





1 



V'/lo'll 





\flnlr: 





_j 



V/iL-2'11-1 







hl/^Rom 

h\'^Ru„ 
h},'^R2,„ 



1/2 



(C.l) 



This new symmetric matrix is called A(t), and R(t) is the matrix of its eigenvectors. The A(t) matrix has the same eigenvalues 
k^{T) as A(t), although the eigenvectors R(t) and R(t) are different but related by R(r) - G(t) R(t). These symmetric matrixes 
are easier to diagonalize numerically. Eigenvectors are computed by the recurrence relation: 

Romir) = 1 

Rl,„(T)^{l-(^(T))lk,„{T) 



1 



RlmiT) = — — [/z,_l(T)«,_l.„(r) - (/- 1)^„(t)/?/_2,,„(t)] 
lk,n{T) 

where, compared to Roberge (1983), w(t) is a r-dependent effective albedo including line absorption. 



(C.2) 



Appendix D: Inverse and derivative of R(t) 

Here we show how R^H''") is computed. Unfortunately, A(t) is not a symmetric matrix, so that R"'(''") 9^ However, we 

can apply the same method as above to turn R"'(t) into R^(t). Since A(t) is symmetric, the matrix formed with its eigenvectors 
is orthogonal. Thus, using the same notations, we have: 



K\t) R(t) = J(t) = (G(t) K{T)f G(t) R(r) 



(D.l) 



where J(r) is a diagonal matrix with 7//(t) 
R-'(r) = r'(T)R^(T)G2(T) 



R;(t) elements. Hence: 



(D.2) 



The inclusion of the depth dependence in the spherical harmonics method unfortunately forces to calculate the R(t) derivative 
respect to t. Ideally, we could start to derivate the Rimij) recurrence relations shown in Eq ( IC.2l i to get: 



Left-multiplication by a diagonal matrix multiplies rows by a constant, and right multiplication multiplies columns. 



= 914.26 A-l = 



10"' 



10" 



10"' 
10-^' 

10-'' 
JO-" 

10-'' 

10-^' 

10"^' 
10-^' 
10-^' 



10 



Af 

f 

-3 



10-' 



10- 




Fig.B.l. Comparison of f and Af for l-Q and A - 914.26 A. The abscissa corresponds to t,, for the upper scale and to T/,„f for 
the lower scale. 
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Im 



io'k„, + (1 - a))k'^ 

km 



k' 1 



with 

h',{T) = -(2Z + 1) W{T)(r,{T) + w(T)tr;(T)] 



(D.3) 



(D.4) 



Unfortunately, w', crj and A:^, have to be computed also numerically, which is quite unstable in the most external cloud positions 
due to the large variations of at line wing wavelengths (where the line opacity becomes comparable to the dust opacity) 
compared to deeper inside the cloud where at the same wavelength becomes saturated (the dust opacity becomes insignificant 
respect to the line opacity). Besides, a symmetric difference scheme does not provide satisfactory results because only 
gives an approximation to o)' ?A t - which, in general, is not t„. We solved this problem by derivating directly the 

computed values of R{t). To avoid irregular steps in r, a second order polynomial was fit to RimiTi-i), Rimini) and Rim(Ti+2), and 
the value of its analytical derivative was then used. The resulting derivative R'(t) is smooth enough to be applied in the numerical 
computation. 

Appendix E: Mean radiation field intensity 



In section [341 we deduced the simple form that the mean intensity takes in the spherical harmonics method, i.e. /(t) - /o(t). 
However, in some cases of astrophysical interest (e.g. a two sides asymmetrically illuminated cloud) one needs to distinguish the 
fraction of radiation field coming from each side of the cloud. In this case, two half sums have to be computed. Here we give the 
analytical expressions that we use to compute /"^(r). For radiation coming from the t = side we have: 



1 fo 1 ^ r" 

And for radiation coming from the t - T,„ax side we have: 
^^(r) = ^^ I{T,iJ.)dfi^^J]i2l+l)fiiT) P,iti)dti 

r+l 

If we define Qi - Pi(m) djj, with: 



(E.l) 



(E.2) 



n 1^0 

Qi - } I even and > 

V /+! 



(E.3) 



I odd 



parity gives Piijj) dfd - (-1)' Qi - -Qi (using Qi -0 for I even). Inserting Eq. (|29T l in Eqs ( IE. II ) and ( IE. 21 ) we get: 
r(T) = 2 Z '^"'^^^ " ^'-^^'^^ ^ ' + ^)QlRhn{T) + - ^ {C„ EliT) + DliT)) 1 - ^(2/ + l)Q,Ri,n(r) 

m=-M V ;=1 J m=\ V /=1 ; 

l ' ( ^ ^l'*^ ( ^ 

J\r) = - 2 (C,„ £,;(t) - D„(T)) 1 + ^(2/ + 1) a /?,,„,(T) + - ^ (C„, ^^(t) + D^t)) 1 + ^(2/ +\)Q, RUt) 



2 . 

m=-M V l=\ ) m=\ V l=\ 

Taking into account the fact that <3/ = if Hs even, and Ri-„, - -Rim if / is odd, we now define (for m > 0) 
Sm{T)^Y,^2l+l)Q,Ri,„{T) 



I odd 



to write: 



-1 M 

r(T) = - ^ (c,„ £,;(t) - djt)) (1 + s„(t)) + - ^ (c„ £,;(t) + d:„(t)) (1 - S,„{t)) 

m=-M in=l 
-1 M 

riT) = - (C,„ E-(t) - D~{t)) (1 - S„{t)) + - 2 (C„ £:(t) + d:„(t)) (1 + S„{t)) 



(E.4) 



(E.5) 



(E.6) 



(E.7) 



(E.8) 



m=-M 



m=I 



Therefore, the fraction of the mean intensity coming from each side of the cloud can be easily determined at each depth. The 
resulting 7^(t) values can then be used to evaluate the escape probably of any FUV photon emitted within the cloud. 



